///\todo upgrade

//////////////////////////////////////////////////////////////////////////////////
//                                                                              //
//  This file is part of the buola project (https://code.google.com/p/buola/).  //
//                                                                              //
//  Copyright(c) 2007-2012 Xavi Gratal                                          //
//  gratal AT gmail DOT com                                                     //
//                                                                              //
//  Buola is free software: you can redistribute it and/or modify               //
//  it under the terms of the GNU General Public License as published by        //
//  the Free Software Foundation, either version 3 of the License, or           //
//  (at your option) any later version.                                         //
//                                                                              //
//  Buola is distributed in the hope that it will be useful,                    //
//  but WITHOUT ANY WARRANTY; without even the implied warranty of              //
//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the               //
//  GNU General Public License for more details.                                //
//                                                                              //
//  You should have received a copy of the GNU General Public License           //
//  along with buola.  If not, see <http://www.gnu.org/licenses/>.              //
//                                                                              //
//////////////////////////////////////////////////////////////////////////////////

#ifndef _BUOLA_MAT_CGAUSSIANPROCESS_H_
#define _BUOLA_MAT_CGAUSSIANPROCESS_H_

#include <buola/mat.h>
#include <buola/mat/cholesky.h>
#include <buola/mat/prettyprint.h>

namespace buola { namespace mat {

template<typename tScalar>
tScalar cov_fun(tScalar pSigma,tScalar pFac,tScalar pDist2)
{
    return pSigma*exp(-pFac*pDist2*0.5);
}

template<typename tScalar>
class CGaussianProcess
{
public:
    CGaussianProcess()
    {}
    
    template<typename tPoints,typename tValues>
    CGaussianProcess(tPoints &&pPoints,tValues &&pValues,tScalar pFac=1.0,tScalar pNoiseVar=1.0,tScalar pOutputVar=1.0)
    {
        Train(std::forward<tPoints>(pPoints),std::forward<tValues>(pValues),pFac,pNoiseVar,pOutputVar);
    }

    template<typename tPoints,typename tValues>
    void Train(tPoints &&pPoints,tValues &&pValues,tScalar pFac=1.0,tScalar pNoiseVar=1.0,tScalar pOutputVar=1.0)
    {
        mFac=pFac;
        mNoiseVar=pNoiseVar;
        mOutputVar=pOutputVar;
        mPoints=std::forward<tPoints>(pPoints);
        mValues=std::forward<tValues>(pValues);
        mMean=mean(mValues);
        for(int i=0;i<mValues.size();i++) mValues[i]-=mMean;
        
        mCov.Resize(mPoints.Rows(),mPoints.Rows());

        for(int i=0;i<mPoints.Rows();i++)
        {
            for(int j=0;j<mPoints.Rows();j++)
            {
                double lDist2=distance2(mPoints(i,nAll),mPoints(j,nAll));

                mCov(i,j)=cov_fun(mOutputVar,mFac,lDist2);
                if(i==j) 
                    mCov(i,j)+=mNoiseVar;
            }
        }
        
        mLLT=cholesky(mCov);
        mKernel=mLLT.Solve(mValues);
    }

    template<typename tPoint>
    tScalar F(const tPoint &pPoint)
    {
        CRow<tScalar> lCov(mPoints.Rows());

        for(int i=0;i<mPoints.Rows();i++)
            lCov[i]=cov_fun(mOutputVar,mFac,distance2(mPoints(i,nAll).T(),pPoint));

        return lCov*mKernel+mMean;
    }

    template<typename tPoint>
    tScalar V(const tPoint &pPoint)
    {
        CVec<tScalar> lCov(mPoints.Rows());

        for(int i=0;i<mPoints.Rows();i++)
            lCov[i]=cov_fun(mOutputVar,mFac,distance2(mPoints(i,nAll).T(),pPoint));

        return mOutputVar-lCov.T()*mLLT.Solve(lCov);
    }
    
    tScalar LogLikelihood()
    {
        return -0.5*mValues.T()*mKernel-sum(log(diag(mLLT.L())))-0.5*mCov.Rows()*log(2*M_PI);
    }

    CMat<tScalar> mCov;
    CVec<tScalar> mKernel;
    CVec<tScalar> mValues;
    CMat<tScalar> mPoints;
    tScalar mFac;
    tScalar mNoiseVar;
    tScalar mOutputVar;
    tScalar mMean;
    CLLT<tScalar> mLLT;
};

/*
mat::CVec gaussian(const mat::CVec &pMu, const mat::CMat &pSigma)
{
    mat::CRowVec lRet=mat::random_normal(1,pMu.Length());
    
//    lRet=lRet*chol(pSigma+mat::identity(pMu.Length(),pMu.Length())*0.00000001);
    
    return lRet.T()+pMu;
}
*/


/*namespace mat*/ } /*namespace buola*/ }

#endif
